// FE_solver_2D_Elasticity.cpp : This file contains the 'main' function. Program execution begins and ends there.
//

#include <iostream>
#include "FEM2DVector.h"
#include "TimeDiscrete.h"
#define PI 3.1415926535897
double lambda(CPoint p,double t)
{
	return 1.0;
}
double mu(CPoint, double t)
{
	return 2.0;
}
double f1(CPoint p, double t)
{
	return -(lambda(p, t) + 2 * mu(p, t))*(-pow(PI, 2)*sin(PI*p.getX())*sin(PI*p.getY()))
		- (lambda(p, t) + mu(p, t))*((2 * p.getX() - 1)*(2 * p.getY() - 1)) - mu(p, t)*(-pow(PI, 2)*sin(PI*p.getX())*sin(PI*p.getY()));
}
double f2(CPoint p, double t)
{
	return -(lambda(p, t) + 2 * mu(p, t))*(2 * p.getX()*(p.getX() - 1)) -
		(lambda(p, t) + mu(p, t))*(pow(PI, 2)*cos(PI*p.getX())*cos(PI*p.getY())) - mu(p, t)*(2 * p.getY()*(p.getY() - 1));
}
double u1(CPoint p, double t)
{
	return 0;
}

double u2(CPoint p, double t)
{
	return 0;
}
double coe_fun(CPoint p, double t)
{
	return 1;
}


int main()
{
	int unknown_size = 2;
	CFEM2DVector fem(unknown_size);	

	//read mesh from file
	//the mesh file contains the P,T,Pb,Tb,boundary nodes,boundary edges, and the type of the element
	fem.readMesh("..\\mesh\\cha5_exa1_linear.neu");

	//The initializing process need mesh's information, so the function readMesh should be called first
	fem.Init();

	//the basis type should be consistent with the mesh file
	BASIS_TYPE basis_trial_type = TWO_DIM_LINEAR;
	BASIS_TYPE basis_test_type = basis_trial_type;
	fem.assignBoundNodesValue(0, 0, Dirichlet, u1);
	fem.assignBoundNodesValue(0, 1, Dirichlet, u2);


	SpMat A1 = fem.assembleMatrix2D(lambda, 0, basis_trial_type, 1, 0, basis_test_type, 1, 0);	
	SpMat A2 = fem.assembleMatrix2D(mu,     0, basis_trial_type, 1, 0, basis_test_type, 1, 0);
	SpMat A3 = fem.assembleMatrix2D(mu,     0, basis_trial_type, 0, 1, basis_test_type, 0, 1);
	SpMat A4 = fem.assembleMatrix2D(lambda, 0, basis_trial_type, 0, 1, basis_test_type, 1, 0);
	SpMat A5 = fem.assembleMatrix2D(mu,     0, basis_trial_type, 1, 0, basis_test_type, 0, 1);
	SpMat A6 = fem.assembleMatrix2D(lambda, 0, basis_trial_type, 1, 0, basis_test_type, 0, 1);
	SpMat A7 = fem.assembleMatrix2D(mu,     0, basis_trial_type, 0, 1, basis_test_type, 1, 0);
	SpMat A8 = fem.assembleMatrix2D(lambda, 0, basis_trial_type, 0, 1, basis_test_type, 0, 1);
	

	//SpMat A = fem.matrix_hstack(fem.matrix_vstack(A1 + 2 * A2 + A3, A6 + A7),fem.matrix_vstack(A4+A5,A8+2*A3+A2));
	SpMat A = fem.matrix_vstack(fem.matrix_hstack(A1 + 2 * A2 + A3, A4 + A5), fem.matrix_hstack(A6 + A7, A8 + 2 * A3 + A2));

	Rsv b1 = fem.assembleVector2D(f1, 0, basis_test_type, 0, 0);
	Rsv b2 = fem.assembleVector2D(f2, 0, basis_test_type, 0, 0);
	

	Rsv b=fem.vector_stack(b1,b2);	
	fem.ExportMat("A_before.txt", A);

	fem.treatBoundaryConditions(coe_fun, 0, A, b);
	fem.ExportMat("A_after.txt", A);


	Rsv x = fem.solveLinearEqua(ITERA_BiCGSTAB, A, b);

	fem.Export2VTK("result.vtk",x);
	


}

